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Abstract.  The  phase  transition  from  liquid  to  ordered  state  has  been  explored  in  Molecular 
Dynamics  simulations  of  confined  one-component  plasma  clouds  (as  in  ion  traps).  The  latent 
heat  and  specific  heat  have  been  extracted  and  compared  with  Coulombic  matter.  The  transition 
temperature  is  found  to  be  lower  for  finite  systems  than  for  infinite  matter,  and  the  difference 
depends  on  the  number  of  confined  particles. 


A one-component  plasma  (OOP)  is  known  to  have  a transition  from  a liquid  to  a 
solid  phase  at  a temperature  corresponding  to  T = 170-180  [1],  where  F = (qVa^J/kT, 
and  au,?  is  the  Wigner-Seitz  radius  of  the  volume  occupied  by  one  ion.  The  behavior 
of  the  latent  heat  is  discussed  in  the  literature,  (though  the  one  plot  of  specific  heat 
that  has  been  published  [2]  omits  the  peak  in  the  specific  heat  for  the  transition). 

The  configuration  of  finite,  confined  piasmas  (fewer  than  100,000  ions)  shows  a 
different  type  of  ordering  from  that  in  infinite  matter  at  low  temperatures:  concentric, 
equidistant  shells  with  2-dimensional  ordering  within  the  shells,  instead  of  a bcc 
structure  [3],  The  question  of  how  sharp  these  transitions  are  in  the  finite  systems, 
and  where  they  occur,  has  not  previously  been  explored. 

Molecular  Dynamics  simulations  were  carried  out,  first  for  a system  of  10,000  ions 
confined  in  a static,  isotropic  harmonic  well,  with  a confining  force  of  F = -Kr.  The 
system  was  cooled  gradually  to  a temperature  that  corresponds  to  T = 30,000.  The 
outer  two  shells  for  this  system  are  shown  in  Fig.  1 together  with  the  shell  structure. 

The  potential  energy  of  the  system  was  extracted  after  equilibrium  had  been 
established.  The  potential  energy  is  given  by 

U = l/2K£r2+Zq2/rij  (1) 

/ /./>»' 


where  r.  is  the  radial  distance  from  the  origin  and  r.  is  the  distance  between  two  ions 
and  q is  the  charge  of  the  ions. 

Then,  temperature  increases  were  applied  by  two  different  methods.  In  one,  the 
temperature  was  simply  increased  by  scaling  up  the  velocities,  and  the  system  was 
allowed  to  equilibrate  for  some  40  plasma  periods  (with  about  120  time  steps  in  the 
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FIGURE  1.  The  upper  pari  of  the  figure  shows  the  outer  layers  of  10,000  ions.  On  the  left  is  a very 
cold  system  with  the  outer  shell  shown  in  black,  and  the  second  shell  in  grey.  On  the  right  is  a system 
in  the  liquid  regime  with  no  structure,  with  the  selection  of  particles  in  the  same  radial  regions.  The 
lower  part  of  the  figure  shows  the  radial  density  for  three  values  of  the  temperature.  The  solid  curve  is 
a frozen,  well-ordered  system,  the  large  points  arc  in  a region  where  the  temperature  is  somewhat 
above  the  melting  point  and  the  small  points  arc  in  the  liquid  regime. 

simulation  per  plasma  period).  Then  the  following  -100  plasma  periods  were  used  to 
establish  the  new  partition  between  kinetic  energy  (temperature)  and  potential  energy. 
Alternatively,  the  temperature  was  increased  in  steps  by  assigning  a random  velocity 
to  each  particle  picked  from  a Boltzmann  distribution  with  the  appropriate  mean 
temperature,  and  then  allowing  the  system  to  settle.  This  latter  method  seemed  to 
allow  the  use  of  somewhat  shorter  time  periods  to  establish  a new  equilibrium.  The 


two  procedures  gave  very  similar  results.  These  times  for  equilibration  and  for 
establishing  the  averages  were  arrived  at  after  some  experimentation  — a shorter  time 
gave  non-reproducible  results.  Attempts  were  also  made  at  similar  calculations  in 
which  the  systems  was  cooled,  rather  than  heated.  Just  as  in  laboratory  experiments  to 
determine  the  melting  or  freezing  temperatures  of  materials,  only  the  melting  process 
seemed  to  give  consistently  reproducible  results. 

The  results  are  shown  in  Fig.  2.  At  low  temperatures,  since  the  particles  are 
oscillating  about  their  equilibrium  positions,  the  energy  increments  are  distributed 
equally  between  kinetic  and  total  energies.  At  a temperature,  corresponding  roughly 
to  the  above  mentioned  melting  point  in  the  infinite  OCP,  there  is  a step  in  potential 
energy,  corresponding  to  the  latent  heat  of  melting.  These  data  are  displayed  in  a 
slightly  different  form  in  Fig.  3 where  the  upper  part  of  the  figure  shows  the  potential 
energy  divided  by  the  kinetic  energy  for  both  the  10,000  ion  simulation  and  the 
infinite  OCP,  the  latter  from  [4].  The  lower  part  of  the  figure  shows  the  specific  heat, 
taken  as  the  increments  in  temperature  divided  by  the  corresponding  increments  in 
total  energy.  There  appears  to  be  a small  difference  in  the  transition  energy,  the  point 
where  the  specific  heat  has  its  maximum  value,  in  that  the  peak  in  the  10,000  ion 
system  occurs  at  -20%  lower  temperature:  T = 210,  instead  of  173. 


FIGURE  2.  The  potential  energy  in  a 10,000  ion  system  is  plotted  as  squares  against  the  kinetic 
energy,  in  the  same  units  of  q2/a(n.  At  low  temperature,  the  zero  in  potential  energy  is  chosen 
arbitrarily  as  the  value  at  T = 0.  For  T < 0.005  the  particles  oscillate  about  their  equilibrium  positions 
and  the  total  energy  is  divided  equally  between  the  kinetic  and  potential  energies.  The  deviation  in  the 
potential  energy  at  T = 0.005  corresponds  to  the  transition  from  an  ordered  to  a disordered  structure. 
At  the  bottom,  the  difference  between  the  potential  and  kinetic  energies  is  plotted. 
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FIGURE  3.  The  upper  pari  ol  ihe  Figure  shows  the  total  energy,  U',  divided  by  the  temperature  for 
infinite  Coulombic  matter  [4]  as  circular  dots,  compared  to  the  simulations  for  10,000  ions  as  triangles. 
The  origins  on  the  vertical  scale  arc  displaced.  On  the  right  (low-temperature)  part  of  the  curve  the 
ratio  is  constant  lor  both,  as  expected.  The  lower  part  of  the  figure  shows  the  specific  heat  C derived 
from  the  same  sets  of  data. 

Next,  the  different  measures  of  ordering  were  investigated: 

• the  width  of  the  shells, 

• the  related  maxima  in  the  radial  densities  for  each  shell, 

• the  heights  of  the  maxima  in  the  2-dimensional  correlation  function  g(r) 
within  the  shells, 

• the  heights  of  the  maxima  in  the  3-dimensional  g(r)  over  the  entire  volume  of 
the  cloud. 
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All  of  these  showed  the  maximum  rate  of  change  near  the  transition  temperature  in 
the  specific  heat,  and  at  a temperature  lower  than  the  transition  temperature  for  the 
infinite  OCP.  This  consistency  is  not  surprising,  since  the  ordering  is  the  source  of 
the  latent  heat.  An  example  is  shown  in  Fig.  4. 


FIGURE  4.  In  (a)  the  3-dimensional  correlation  function  g(r)  is  plotted  for  three  values  of  T,  in  (b) 
the  ratio  of  the  maximum  in  g(r)  to  the  following  minimum  is  shown,  also  as  a function  of  T,  while  in 
(c)  the  derivative  of  this  ratio  is  plotted. 
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In  order  to  investigate  the  question  of  whether  the  lowering  of  the  transition 
temperature  is  real  or  not,  similar  calculations  were  carried  out  for  1000  and  100  ions. 
These  are  displayed  in  Fig.  5.  It  is  clear  that  the  peaks  in  specific  heat  are  broader 
and  weaker  in  the  smaller  systems.  For  both  1000  and  100  ions  the  calculations  were 
repeated  several  times  in  order  to  ascertain  that  the  effects  were  reproducible.  As  is 
shown  in  the  lower  part  of  Fig.  5,  there  indeed  appears  to  be  a systematic  trend  in  the 
transition  temperature  becoming  lower  for  the  smaller  clusters. 

Such  a trend  appears  to  have  a simple  explanation,  in  that  ions  on  the  surface  do  not 
see  nearest  neighbors  on  the  outside,  and  thus  their  mobility  is  not  constrained  in  the 
same  way  as  for  ions  in  the  interior.  Such  effects  are  known  to  occur  in  atomic  and 
molecular  clusters  (for  instance  for  gold  atoms  [5]),  but  there  the  interactions  are 
short  range.  One  might  expect,  that  with  ionic  clusters  where  the  interaction  is  purely 
through  the  long-range  Coulomb  force,  the  effect  might  have  a different  character. 
The  lower  part  of  Fig.  5 shows  that  the  transition  temperature,  as  a function  of  the 
fraction  of  ions  that  are  in  the  surface  layer,  seems  to  obey  a simple  dependence. 

For  real  ion  traps,  the  geometry  considered  here  is  close  to  that  in  Penning  traps, 
where  for  the  rotating  cloud  the  effects  of  rotation  at  low  temperature  may  be 
transformed  away  by  the  choice  of  the  coordinate  frame.  The  above  simulations  may 
therefore  be  expected  to  be  directly  tested  there.  The  question  of  the  transition 
temperatures  in  ions  confined  in  rf-traps  (Paul  traps,  and  other  configurations)  is  not 
necessarily  so  clear.  The  fact  that  the  imposed  rf  motion  does  not  cause  appreciable 
heating  of  a cold  plasma  has  been  shown  [6]  but  no  further  studies  have  been  made. 
Therefore,  a simulation  was  carried  out  for  1000  ions  in  rf  confinement  as  described 
in  [6],  In  these  systems  a quasi -temperature  is  defined,  that  includes  only  motion  that 
is  aperiodic  in  the  rf  field:  only  motion  that  corresponds  to  displacements  in  complete 
rf  periods  are  considered.  The  kinetic  energy  in  the  rf  motion  can  be  orders  of 
magnitude  higher.  For  calculations  of  specific  heat,  the  potential  energy  of  the 
system  also  needs  to  be  known,  and  this  clearly  oscillates  during  the  rf  motion.  Two 
methods  were  tried: 

a)  averaging  the  potential  energy  over  complete  rf  cycles,  or 

b)  averaging  the  coordinates  of  each  ion  and  computing  the  potential  energy 
corresponding  to  these  average  coordinates. 

The  two  methods  gave  essentially  identical  results.  One  might  have  expected  that 
because  the  method  of  confinement  involves  a continuous  driven  oscillation  of  all 
particles,  plasmas  confined  in  such  a field  might  perhaps  be  more  fragile  and  melt 
more  easily  (at  lower  temperatures)  than  those  confined  is  static  field.  However,  the 
specific  heat  from  the  simulations  in  the  rf  field  peaked  at  a value  of  F that  was 
consistent  with  the  static  calculations  and  no  such  effect  was  observed. 

The  values  of  specific  heat  much  above  the  transition  temperature  could  not  be 
extracted  from  the  rf-confined  system  because,  as  was  shown  in  [6],  at  these 
temperatures  there  is  increasing  self-heating,  "rf  heating",  from  the  coupling  of  the  rf 
motion  into  the  random  heat  of  the  plasma  and  to  obtain  meaningful  results  for  these 
higher  temperatures  will  require  further  study.  However,  it  is  clear  that  the  rf  motion 
does  not  appear  to  make  systems  confined  in  such  an  environment  any  more  fragile 
with  respect  to  melting  than  those  confined  in  static  fields. 
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FIGURE  5.  The  upper  part  of  the  figure  shows  specific  heats  from  simulations  of  10,000,  1000  and 
100  ions  (circles,  squares,  and  stars).  The  lower  part  of  the  figure  shows  the  transitions  temperatures 
(temperatures  at  which  there  is  a peak  in  the  specific  heat)  as  a function  of  the  number  of  ions  in  the 
surface  layer  for  these  systems.  The  line  corresponds  to  T = T,)(..„  x (I  - 0.96  x fraction). 
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